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Abstract 

The tidal disruption of a star by a massive black hole (MBH) is thought to produce a transient 
luminous event. Such tidal disruption events (TDEs) may play an important role in the detection 
and characterization of MBHs and probe the properties and dynamics of their nuclear stellar clusters 
(NSCs) hosts. Previous studies estimated the recent rates of TDEs in the local universe. However, the 
long-term evolution of the rates throughout the history of the universe has been hardly explored. Here 
we consider the TDE history, using evolutionary models for the evolution of galactic nuclei. We use a 
ID Fokker-Planck approach to explore the evolution of MBH-hosting NSCs, and obtain the disruption 
rates of stars during their evolution. We complement these with an analysis of TDEs history based on 
N-body simulation data, and find them to be comparable. We consider NSCs that are built-up from 
close-in star-formation (SF) or from SF/clusters-dispersal far-out, a few pc from the MBH. We also 
explore cases where primordial NSCs exist and later evolve through additional star-formation/cluster- 
dispersal processes. We study the dependence of the TDE history on the type of galaxy, as well as 
the dependence on the MBH mass. These provide several scenarios, with a continuous increase of the 
TDE rates over time for cases of far-out SF and a more complex behavior for the close-in SF cases. 
Finally, we integrate the TDE histories of the various scenarios to provide a total TDE history of the 
universe, which can be potentially probed with future large surveys (e.g. LSST). 


1. INTRODUCTION 

A tidal disruption event (TDE) by a massive black hole 
(MBH) occurs when a star comes closer than its tidal ra¬ 
dius (approximately given by rt = (M,/m*)^/^r*; where 
M, is the mass of the MBH, m* is the mass of the star 
and r* is its radius). The star is then pulled apart by the 
tidal forces by the MBH, part of its materi al is ejecte d 
away and part is accreted on the MBH ('e.g. lRee(^ri988l l. 
The explosive disruption itself, and the later accretion on 
the MBH likely give rise to transient energetic phenom¬ 
ena potentially observable (with more than ten potential 
detec tions to date ) with current and future telescopes 
(e.g. iCezaril[2014( 1. 

MBHs are thought to exist in a significant fraction of 
all galactic nuclei. Due to their high mass (10^ —IO^Mq), 
they dominate the gravitational potential in the nuclear 
star clusters (NSCs) which typically host them (at least 
at the 10® — 10® Mq MBH mass regime); the orbits 
of stars in the central parsecs around MBHs are there¬ 
fore governed by the potential of the MBHs. The long¬ 
term dynamical evolution of the stars in the region of 
influence of the MBH is also determined by their inter¬ 
actions and mutual gravitational scatterings with other 
stars in the NSC. These can be modeled using a Fokker- 
Planck (FP) approa ch to describe the diffu sive behav¬ 
ior of stars in NSC. ((Bahcall and WoIlll976f) . where the 
close approaches of stars to the MBH leading to their 
tidal disruption can be modeled as a sink-term in the FP 
equations. TDE rates calculated following such an ap¬ 
proach typically give rise to rates of the order of 10“® 
up to \0~^yr~^ per galaxy in the local Universe, as- 
suming some stea.dy st ate condition of the NSC (e.g. 
iWang and Merritul2004 and references therein). 

Since the 1980’s tidal disruption of stars by MBHs 
were suggested to occur in galactic nuclei and produce 


lumin ous transient events with unique signatures (jReesI 
1198811 . The signatures can teach us both on the dy¬ 
namics in galactic nuclei as well as on the conditions 
very close to MBHs, and their properties. Such TDEs 
have since been extensively studied theoretically, and 
a significa nt observatio nal effort has been put in their 
detection ((Cezaril I2014D . Many studies have estimated 
the rate s of T DEs l e.g. iMa gorrian and Tremaind 
iMerrittI (|2009ll and iBrockamp et al.l 1201111 in the local 
universe, and several estimates have been derived based 
on direct observations of TDE-candidates in recent years 
(iKhabibullin and Sazono"^ 12014 Ivan Velzen and Farrail 
1201411 . Since the structure and the properties of NSCs 
vary with time, it is expected that the TDE rates will 
change accordingly. Moreover, the change of TDE rates 
can be used as a proxy for characterizing the proper¬ 
ties of NSCs and/or the history of the build-up of MBHs 
and NSCs in the universe. In particular, future deep sur¬ 
veys may allow the detection of TDEs from earlier stages 
of cosmological evolution arising from younger galaxies, 
and potentia lly provide a handle on the TDE history 
(|Gezarill2014fl . Nevertheless, the evolution of TDE rates 
with time and their dependence on the NSC and star- 
formation history and evolution have been little studied. 

In this work, we explore the TDEs history (TDH) by 
studying the rate of tidal disruption of stars in NSCs. 
We simulate the evolution of NSCs in several types of 
galaxies by FP analysis as well as N-body simulations. 
We obtain the TDE rates as a function of their evolu¬ 
tionary history, and find that they are affected by the 
NSC structure and its build-up history. By integrating 
the TDE rates from galaxies of each type we obtain the 
total galaxy averaged TDH of the universe expected to 
be observed. We thereby provide a first global estimate 
of the TDH accounting for NSC build-up, albeit with 
simple models. 
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In the following we begin by a brief description of the 
Fokker-Planck analysis and N-body simulations of evolv¬ 
ing clusters. We describe our approach of evaluating a 
sink term from FP analysis for calculating the tidal dis¬ 
ruption rate. We then present the TDH and its depen¬ 
dence on the MBH mass and the NSC build-up history 
in different types of galaxies. Finally, we consider the 
global, galaxy averaged, TDH of the universe evaluated 
by integrating the TDE rates arising from the different 
types of galaxies (and their correspondent evolutionary 
histories) and relevant MBH masses. We note, that our 
model is based on typical NSCs that undergo basic pro¬ 
cesses and their dynamics are dominated by two-body 
relaxation. The model can be expanded, and future 
works may include other processes and properties in NSC 
such as non uniform stellar mass, self-consistent poten¬ 
tial, non-spherical NSCs etc. Here we provide a first 
study of the history of TDEs using simple models. 

2. METHODS 

In our work we explore the tidal disruption rates in 
NSCs through their evolution in time. We focus on sim¬ 
ple, spherical NSCs that evolve throug h in-situ star for¬ 
mation (SF) (| Aharon and Peretil2015ll . and study other 
cases of NSCs that are built-up through consecutive in¬ 
falls of massive clusters or evolve from a pre-existing stel¬ 
lar cusp. The main analysis method used in our work is 
based on the numerical solution of the FP equation, com¬ 
plemented by N-body simulations of NSCs evolution. 

Our model includes a stellar cluster, the NSC, harbor¬ 
ing a central MBH. It focuses on the stars in the central 
few parsecs of the NSC, and in particular in the range be¬ 
tween the tidal radius - n ^ below which 

regular stars are disrupted by the MBH, and the radius 
of influence where the stellar motions are dominated by 
the MBH potential, defined by 

= GM,/al (1) 

where r*, and M* are the typical radius and mass of stars 
in the NSC, respectively; and ct* is the velocity disper¬ 
sion of stars just outside the NSC. For r > rt the original 
model following iBahcall and Wc^ (jl97(il . hereafter BW) 
assumes the existence of a “thermal bath” which supplies 
stars to the inner region of the galactic nucleus. The stel¬ 
lar orbits within r/j are assumed to be Keplerian in this 
range. The relaxation time, that dominates the timescale 
of the stellar cluster, is defined as 


ln(A) ^ ^ 

Where p is the stellar density, and In (A) is the Coulomb 
logarithm, a factor which is related to the scale of the 
system (ln(A) 10). 

2.1. Fokker-Planck analysis and N-body simulation 

The FP model, extensively used in our work, consists of 
a time and and energy-dependent, angular momentum- 
averaged particle conservation equation. It has the form: 

^i^^ = -AE-i^-FLc{E,T)+FsFiE,T) (3) 


where 

A = ln(A) (4) 

The term F = F[f(F),Ej is related to the stellar flow, 
and plays an important role in the evolution of the stellar 
cluster. It presents the flow of stars in energy space due 
to two-body relaxation, it is defined by: 


FsF(£;,T) = ^(n(E)EoE‘^) ( 6 ) 

at 

n(i?) is a rectangular function, which boundaries cor¬ 
respond to the region where new stars are assumed to 
from; Eq is the source term amplitude; and Fsf is a 
power-law function with a slope a, defining the SF dis¬ 
tribution in phase space. Fo r a detailed descrip t ion o f 
the Fokker-Planck analysis see lAharon and PeretsI (120151 1 
where the same analysis was used for studying the evo¬ 
lution of NSCs with in-situ SF. Note that the sink term 
Flc{E,T) in Eq. [3] is used for evaluating the tidal dis¬ 
ruption rates in our models. This term represents stars 
with energies in the interval {E, E -\- dE) that flow into 
the MBH (mostly due to angular momentum change; 
the sink-term represent the effective loss in energy space ; 

see lBahcall and Wolil977HLightman and Shaoir nil 97711 

and are therefore lost from the system (see lMerrittil2013al 

for detailed overview of loss cone analysis). 

For modeling tidal disruption through N-body sim¬ 
ulations, we follow the same methods and similar as¬ 
sumptions, as well as make the same use of the 
sam e code in iP erets and Mastr obuono-Battistil (120141 1 
and iMastrobuono-Battisti et al.l (|2014ll . where the TDE 
rates are evaluated. For a detailed method explanation 
along with the initial condit ions of the infalling clusters 
and T DEs estimation, see IMastrobuono-Battisti et Ml 
(|2014l l. In brief, we r un direct N-body s imulations (using 
the i^GRAPE code, iHarfst et al.ll200'i^ of the consecu¬ 
tive infall and merging of a set of 12 single-mass globular 
clusters each with mass of 10® Mq, inspiralling from a 
galactocentric distance of 20pc. 

3. MODELS 

3.1. Formation and evolution of NSCs 

In our study, we explore the TDH through two forma¬ 
tion/evolutionary models of NSCs: in-situ star forma¬ 
tion and cluster infall. The study is done by using the 
FP analysis to simulate an evo lving NSC with an addi- 
tion of extra source term (see lAharon and Per^ 120151 
for the full analysis). In the in-situ SF model we con¬ 
sider two possible cases: NSCs that are entirely built-up 
from the in-situ SF, and NSCs that evolve from an ini¬ 
tial steady state BW distribution (hereinafter “built-up 
SF” origin and “primordial cusp origin”, respectively). 
In both of the scenarios, the initial density profile is de¬ 
termined from an arbitrary initial distribution function 
(DF). For the bound stars inside the radius of influence 
in the primordial cusp scenario the DF takes the form of 


F = J dE lf{E, t)—^, - f(E , t) j (max(E, E )) 

( 5 ) 

The addition of the source term that represents the SF 
has the form: 


tolu 
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f{E,to) oc (corresponding to the BW steady-state 
cusp which has the form of n oc The normaliza¬ 

tion of the number density at the radius of influence for 
each scenario is given in table 13.11 We consider different 
masses for the hosted MBH in their centers, where the 
density profiles of the NSCs are n ormalized in the same 
way as done in I Wang and Merritt! (2004), where the pa¬ 
rameters were taken from typical galaxies (e.g. NGC 
4551/4621). The velocity dispersion (which is taken into 
account in the calculation of the radius of influence) 
is evaluated from the M-Si gma relation: M, oc 
(|Merritt and Ferraresell200lh . 


1 MBH mass [M©] 

1 5 X 10*= 

1 1 X 10“ 

1 5 X 10“ 

1 1 X 10' 1 

1 density profile at 

1 6 X 10"* 

1 1 X 10"* 

1 4 X 10"* 

1 1 X 10“ 1 


TABLE 1 

Density profile at of NSCs for different MBH masses 


In the built-up SF scenario the initial DF is evaluated 
according to the term Esf (Eq. [HI), where the initial 
number density is equal to the added number density 
every time step. For both of the scenarios. The chosen 
slope of th e SF function i s motivated by the observed 
power-law (|Do et al.l 120091 : iBartko et al.l feoOQD distribu¬ 
tion of young stars observed in the young stellar disk in 
the GC. The DF of the unbound stars (—oo < r < r^) 
has the form of a Maxwellian distribution. We study the 
TDE rates throughout the evolution of NSCs and For 
each scenario we tested several types of SF history cor¬ 
responding to the galaxy type studied (see Section 13.21) . 
Though we use the FP approach for our modeling of the 
TDH, we also complement our analysis with N-body sim¬ 
ulations of the cluster-infall scenario. As discussed later 
on, we find that models of in-situ SF in the outer regions 
of an NSC evolve very similarly to those of cluster infall, 
since the clusters shed their material in the same regions, 
effectively providing a source term similar to that aris¬ 
ing from SF. Therefore our models of in-situ SF in the 
outer regions of an NSC effectively well capture models 
of cluster-infall, at least in terms of the TDE rates under 
the assumptions made in our study. 

3.2. Star formation history in different type of galaxies 

We followed the evolution of several types of galax¬ 
ies with different SF scenarios. For each type of galaxy 
and SF scenario we evaluated the TDE rates and their 
evolution. We considered the following simplified cases: 
(1) An elliptical galaxy with SF occurring only at the 
early stages of the evolution of the NSC; and (2) spi¬ 
ral galaxy with either continuous SF or repeating bursts 
of SF occurring throughout the evolution of the NSC to 
present days. The correspondenc e between galaxy type 
and it s typical SF history follows iMadau and DickinsonI 
(|2014ll . We used a simplified model to describe the non- 
continuous SF star-burst scenario, where we assume the 
bursts occur at equally spaced time intervals (Gyr), at 
ten times the rate of the continuous SF models, but only 
for 100 Myrs (as to have the same total SF averaged 
rate). We summarize the types of galaxies and the SF 
scenarios simulated in our work in Table [D Note, that 
we also study different scenarios for the regions where 
SF may occur; either in the range between 2.0 — 3.5pc 


or in the range 0.05 — O.lpc from the MBH. The ranges 
that were chosen for SF in our work are motivated by the 
observed existence of a possibly star-forming circumnu- 
clear gaseous disk in the range of 2-5pc in the Galactic 
center (similar young stel lar populations a re observed in 
other galactic nuclei, e.g. iSeth et ;ini2006f) . Most of the 
presented SF scenarios are based on these observations. 
The close-in SF region (< O.lpc) is motivated by the ex¬ 
istence of a nuclear stellar disk of young stars very close 
to the MBH in the Galactic center (e.g. iDo et all 120091 : 
IBartko et al.ll2009ll . The averaged rates of SF were set 
to 10“'^ — 5 X 10“^yr“^, motivated by the observed num¬ 
ber of stars in the central parsecs of galactic nuclei, as to 
obtain a total mass of ~ 10® — 5 x IO’^Mq comparable to 
the masses of NSCs deri ved from observat ions of the GC 
and extragalactic NSCs (|Seth et al.ll2006l) . 

In addition, we note that the range of of 2.0 — 3.5pc is 
also representative of the range where a disrupted cluster 
dispenses its stars close to the M BH, thereby building- 
up the NSC stellar population (e.g. lAntonini et al.l[2012l : 

iPerets and Mastrobuono-Battistill20l4 the exact cluster 

tidal radius depends on the MBH and NSC mass and 
could have a slightly larger range). Finally, we also con¬ 
sidered different rates of SF in the NSC to access their 
overall effects on the TDH. All models used are summa¬ 
rized in Table 1. Note that bursting SF models where 
SF occurs in the outer regions also serve as good proxies 
for the cluster infall model, as we discuss later on. 

We simulated the evolution of each NSC (described in 
Table (H) with its corresponding TDE for four different 
masses of MBHs: 5 x 10 ®Mq, 1 x 10 ®Mq, 5 x 10 ®Mq 
and 1 X IO^Mq in order to evaluate a general form of 
TDEs dependency on MBH mass. For each studied mass 
of MBH the appropriate relaxation time corresponding 
to the velocity dispersion (evaluat ed from M-Sigma re¬ 
lation, (MerHtE^nT^rraresillQOll) was determined, and 
used to produce the primordial clusters considered in the 
relevant models. The rates of tidal disruptions were then 
evaluated throughout the NSC evolution. The evolution 
of NSCs hosting MBHs of higher masses is not discussed 
in this work, since the tidal radius for MS stars falls 
below the Schwarzschild radius {2GM,jc^). Therefore, 
MBH with higher masses would not contribute to the 
TDE rates of typical main-sequence stars. Lower mass 
MBHs may also exist, though there is little observational 
evidence for their current existence, but this is likely due 
to current observational limitations. Nevertheless, such 
low mass MBHs may give rise to TDEs with significantly 
different observational signatures, and we therefore limit 
our study only to more massive MBHs. 

3.3. Assumptions and Limitations 

The disruption rates of stars by MBHs depend both on 
the NSC properties as well as on several possible phys¬ 
ical processes. These processes include for example the 
prese nce of massive per turbers, that refill the loss cone 
(e.g. IPerets et al.ll2007t) . the evolution of central MBHs 
binaries and the effects of n on spherically symmetric po¬ 
tentials in the nucleus Isee lMerrittJl2013bl for a compre¬ 
hensive review). Some of these processes are efficient 
in increasing the TDE rates, but typically operate only 
on short timescales (e.g. binary MBH mergers), while 
the level of non-sphericity of galactic nuclei is not well 
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Model (Galaxy type) 

SF rate {yr -*-) 

Spatial regions 
of SF (pc) 

Epoch of SF 

NSC Origin 

E (elliptical) 

5 X 10-3 

2.0 - 3.5; 

0.05 - 0.1 

0-1 Gyr 

Primordial BW 
cusp 

0-1 Gyr 

buiid-up 

10-3 

2.0- 3.5 

0-1 Gyr 

Primordial BW 

cusp 

0-1 Gyr 

buiid-up 

(spiral with continuous SF) 

5 X lO-'^ 

2.0- 3.5; 

0.05 - 0.1 

0-10 Gyr 

Primordial BW 

cusp 

0-10 Gyr 

buiid-up 

Sb (spiral with bursts of SF) 

Bursts of 5 X 10-3 

2.0- 3.5; 

0.05 - 0.1 

bursts of 100 
Myr every 1 

Gyr 

Primordiai BW 

cusp 

buiid-up 

Bursts of 10-3 

2.0- 3.5 

Primordiai BW 
cusp 
buiid-up 


TABLE 2 

Models for nuclear stellar clusters and star formation in different types of galaxies. 


known. Our model focuses the TDEs for number of SF 
scenarios including different rates and different ranges. 
A sim ilar TDE history model was suggested bv iMerritl 
(|2009r) whose work was based on different evolution sce¬ 
narios of NSCs which does not include SF. A comparison 
between Merritt ’s work and this study results is discussed 
in section [4. 1.11 

Another important assumption is that the MBHs in 
our models evolve and change their mass. This change 
is done by adding the mass of the disrupted stars to the 
MBH during the evolution of an NSC. We emphasize, 
however, that the growth of the MBH mainly affects the 
TDE rates in the lower rates of SF scenarios and for 
initial low masses of MBH, in which cases the total added 
mass becomes comparable to the initial mass of the MBH 
. In these cases, the TDE rates decline with time (see 
Fig. O, while the TDE rates grow with time for the 
more massive MBHs. 

We continuously change the NSC structure up to the 
radius of influence according to the evolution of the FP 
models. As the mass of the MBH grows we change 
the velocity dispersion cr* assuming the M-cr relation 
holds at any given time. The radius of influence (Eq. 
d]) is then changed according to the updated mass and 
velocity dispersion, and the outer boundary condition 
of the thermal ba t h also change accordingly, following 
IWang and Merritt (l2004f) . We note and caution that this 
is not a trivial assumption; its advantage is in consis¬ 
tently relating the MBH, the velocity dispersion and the 
stellar density outside the NSC, and therefore, by defi¬ 
nition producing a final NSC and environmental config¬ 
uration which is in accordance with the observed M — cr 
relation today (though the relation is defined in the con¬ 
text of the velocity dispersion at larger radii). The dis¬ 
advantage is in using an ad-hoc assumption in which the 
M — crrelation is continuously kept at all times; whether 
this holds is currently unknown. To the best of our 
knowledge previous studie s used either a constant outer 
boundary conditions (e.g . iBahcall and Wolflll976l[l9^ 
iCohn and Kulsrudlll97^ . o r had self-consistent ly evolv¬ 
ing open outer boundaries (iMurphv et ahlflO^ : in this 
latter case the initial conditions determined the velocity 
dispersion and the density in the NSC outer region. Since 
no additional boundary conditions or an outer thermal 
bath were introduced, the conditions in the outer regions 
(where relaxation is slow) hardly changed (see below), 


and the NSCs was effectively described as isolated sys¬ 
tem, unrelated to the larger scale evolution of the galaxy. 
Thereby it could not reproduce a configuration based on 
the M — cr relation today, unless arbitrarily put by hand 
as initial conditions. Given these limitations and the de¬ 
bated understanding of MBH feedback and the origin of 
the M — cr relation, we believe the ad-hoc assumption we 
use is quite reasonable, and compared with other choices 
it introduces a plausible relation between the NSC and 
the larger scale environment consistent with current ob¬ 
servations. 

Our modeling accounts for the NSC stellar population 
(the density profile) up to the radius of influence in the 
evolving NSC, while the stellar population outside the 
region is assumed to be consistent with the large scale 
galactic bulge, as described above. We assume the po¬ 
tential is dominated by the MBH and neglect the contri¬ 
bution from the stellar population of the NSC, i.e. the 
contribution of the stellar component up to the radius of 
influence is neglected, similar to the original BW study. 
Up to the radius of influence this simplifying assumption 
is justified, while far from the MBH it breaks down. Nev¬ 
ertheless, at these regions the relaxation time becomes 
larger than the Hubble time, and the background does 
not change due to the NSC internal evolution. The inter¬ 
mediate regions beyond the radius of influence, in which 
relaxation time can still be significantly shorter than a 
Hubble time are small, and therefore our simplified mod¬ 
eling should still well capture the overall dynamics of the 
NSC. 

For our final integrated TDH results (see Section 031) , 
we extended the range of MBHs masses between 10® and 
5 X 10^ Mq by using the interpolated power-law function 
that describes the dependence between TDE rates and 
MBH mass (Eq. [7|). We do not consider TDEs from 
lower mass MBHs, or TDEs of evolved stars, which could 
be relevant even for MBHs above this mass range. 

Here we only consider simple models of spherical NSCs 
hosting non-evolving MBHs of given masses with relax¬ 
ation dominated by two-body scattering of low mass 
(Sun-like) stars. Other models are beyond the scope of 
this initial study and will be explored elsewhere. 

4. RESULTS & DISCUSSIONS 

4.1. Tidal disruption events history from Fokker-Planck 

analysis 
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4 . 1 . 1 . Cluster infall and star-formation in NSC outskirts 

Stellar clusters infalling into the nuclear regions are 
typically disrupted and shed their stars at scales of a 
few pc from the Galactic center, due to their disrup¬ 
tion by the potential of the MBH and the ex isting NSC 
(e.g. iPerets and Mastrobuono-Battistill2014f) . There is 
also evidence for star-formation, young stellar popula¬ 
tions and dense star-forming clumps in these regions (e.g. 
the nuclear stellar disk in the Galactic center, and young 
stellar populations in many galactic nuclei, iSeth et HI 
[2fM[0^ et al.ll201lll . Together, these suggest that both 
in-situ SF and cluster infall scenarios give rise to a source 
term of stars in regions at a few pc from the MBH. In the 
following we present the results of FP calculations, where 
star-formation/cluster infall scenarios are assumed to in¬ 
troduce new stars, a “source term” in the FP models, in 
the regions between 2-3.5 pc. Note that we use the term 
SF throughout the following discussion, but it effectively 
also represents a source term of stars from infalling clus¬ 
ters. 

As discussed in the previous sections, we have explored 
the TDE rates as a function of the SF rate and charac¬ 
teristics. Figures [Hand [2] show the evolution of the TDE 
rates with time for various models of the NSCs. As ex¬ 
pected, the TDE rates grow with time as the NSC stellar 
population and the density of the NSCs increase due to 
star-formation. The early times TDE rates in the mod¬ 
els of early star-formation (ellipticals) are significantly 
higher than the corresponding rates in the models with 
long term SF (spirals). However, at later times the NSCs 
in both cases saturate to similar levels as the NSCs reach 
similar densities at late times (in cases of the same to¬ 
tal cumulative SF throughout the evolution). The TDE 
rates in NSCs with initial primordial cusps are generally 
higher, since the overall stellar populations are larger. 
The rise in the TDE rates is not linear as in the built- 
up NSCs models as relaxation processes are more pro¬ 
nounced already at early stages. Since typical relaxation 
times are longer than the length of SF bursts in our mod¬ 
els, it is not surprising that the SF bursts models behave 
very similarly to those with continuous SF (with same 
averaged SF rates). Henceforth we omit the results for 
SF burst models, which are essentially the same as those 
in the continuous SF models. 

The high SF during the first Gyr in the elliptical galax¬ 
ies models give rise to a significant increase of the stel¬ 
lar population in the NSC, and correspondingly leads to 
shorter relaxation times, and higher densities. In turn 
the TDE rates rapidly increase during this epoch, due 
to these effects, and the increase becomes shallower only 
at later times after the end of the SF epoch. Neverthe¬ 
less, the rates continue to grow in time even after 10 Gyr 
of evolution as stars from the galaxy slowly diffuse into 
the NSC. In the spiral galaxy models SF is continuous 
throughout the evolution, and the TDE rates similarly 
increase in a more continuous manner. 

In addition to the comparison between TDEs in differ¬ 
ent types of galaxies, we compare the TDE rates for dif¬ 
ferent masses ofMBHs (5x10®M0, 1x10®Mq, 5x10®Mq 
and 1 X 1O^M0). We present the behavior of each sce¬ 
nario and for each of the MBH masses along with the 
evolution of the total mass of the NSC; see Figures O 
and 01 These probe the dependence of the TDE rates 


on the mass of the MBHs. It can be clearly seen that 
the TDE rates decrease with the mas s of the MBH. This 
result has already been discussed by iWang and Merritt 
(|2004|1 . and can be well understood from loss cone anal- 
v sis ('seelMerrittJl2013al l. 

iMerrittI (120091 ) have discussed the evolution of TDE 
rates with time in the case of a non-SF NSC. Though our 
models always include SF, the evolution of models of pre¬ 
existing NSCs in which only an early SF epoch occurs, 
evolve very similarly to the models explored in IMerrittI 
(|2009f) . after the end of the SF. Such models generally 
predict a global decrease of TDE rates over time. . A 
different trend exists in the evolution of TDE rates in 
models with long-term SF, in which the TDE rates slowly 
increase with time as the NSC keeps growing due to SF. 

Fig. 0] shows the TDE rates for the models with lower 
SF rates and low MBH masses. Note that in the models 
with MBHs masses of 5 x 1O®M0 and 10®M0 the TDE 
rates decay with time. In these cases the restive increase 
in the MBH mass significantly increase the relaxation 
time leading to lower TDE rates, even though the total 
masses of the NSCs increase. 

4 . 1 . 2 . Star-formation in the inner regions of NSCs 

In addition to the TDH presented above, correspond¬ 
ing to SF at a distance of 2 — 3.5pc, we also studied 
the TDH in a different scenario, where SF is assumed to 
mainly occur in the inner regions of NSCs. Evidence for 
such more nucleated SF can be found in our own galaxy, 
where a young stellar population is observed in the inner 

O.lpc near the MBH; though it is not clear if such SF 
is typical (such regions can not yet be resolved in other 
galaxies), we consider this possibility for completeness. 
We followed the same conditions used in the previous 
scenarios, but we changed the range of SF to correspond 
to distances of 0.05 — O.lpc. Figure |6| presents the TDEs 
history for SF rate of 5 x 10~^yr~^ and 5 x 
for elliptical and spiral galaxies, respectively. On the 
right panel we notice a different behavior of the TDH 
for elliptical galaxy compared with the outskirt SF sce¬ 
narios: while there is a rapid growth in the disruption 
rates during the epoch of SF, the rates decreases with 
time once the SF process is quenched. Since the stellar 
densities increase on relatively short times as many stars 
are formed close to the MBH, the NSC is initially far 
from its steady state configuration. The inner regions are 
“over-populated” and relaxation processes take time to 
redistribute these stars throughout the NSC. Since stars 
in these regions are more susceptible to being tidally dis¬ 
rupted (larger loss cone), the TDE rates increase rapidly, 
before the NSC structure re-equilibrates. The long-term 
SF occurring in spiral galaxies allows for relaxation pro¬ 
cesses to redistribute the stars, and the TDE rates in¬ 
crease more slowly, as the NSC is built-up. 

4 . 1 . 3 . Global SF in NSCs 

In addition to the TDHs of NSCs evolving through 
an early SF epoch or through continuous SF (de¬ 
scribed in Table ED, we also studied models in which 
the SF follows the global universal SF history (i.e. 
iMadau and Dickinsonll20T^ in such global models early 
SF is more significant, and then it gradually decrease b y 
an order of of magnitude (|Madau and Dickinsonll^l4f) . 
Such models are therefore expected to be more similar to 
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x10“® x10'® 




Fig. 1.— The evolution of TDEs for NSCs with an MBH of 4 x 10 ®Mq for 3 types of galaxies, that are evolved through in-situ SF (left) 
and from a primordial cusp (right). The rates of SF are 10“^ and 10“"^ starsxyear"^ for elliptical and spiral galaxies respectively. The 
total number of stars formed throughout the evolution is the same; NSC with primordial cusps have a larger number of stars as they have 
similar SF histories but they also include the stellar populations of the primordial cusps. 


x10“^ x10“^ 



respectively). 
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Fig. 3.— Evolution of TDE rates (continuous lines) for four different masses of MBHs for model E along with mass of the NSC (dashec 
lines). The cyan vertical dashed line marks the epoch at which the SF stops. NSC masses and the corresponding TDE rates of the sam< 
models are marked with the same colors. Left: TDH for NSC evolved from built-up star formation. Right: TDH for NSC evolved from i 
primordial cusp. 




Fig. 4. — Evolution of TDE (continuous lines) rates for four different masses of MBHs for models S'o (solid lines) and Sb (dashed lines). 
NSC masses and the corresponding TDE rates of the same models are marked with the same colors. The blue dashed lines mark the time 
of the SF bursts (every 1 Gyr). Left: TDH for NSC evolved from built-up star formation. Right: TDH for NSC evolved from a primordial 
cusp. 


the early SF epoch models (“elliptical galaxies”). In Fig. 
[3 we show the results from thee models, which confirm 
these expectations; these models show similar trends as 
those in Fig. [71 

4.2. Tidal disruption events history using N-body 
simulation data 

Following the cluster infall study (described in 
iPerets and Mastrobuono-Battistil 1201411 . we evaluated 
the disruption rates of stars based on data obtained from 
N-body simulations. This provides us with a complimen¬ 
tary approach to the FP calculations, through a differ¬ 
ent type of analysis. As we show in the following the 
FP calculation for models with SF in the NSC outskirts 
and the calculations based on N-body data are gener¬ 
ally comparable. Note that the N-body calculations are 
computationally expensive, and therefore we show only 
a single case, while our general main results (shown in 


the previous sections) are based on the FP calculations. 
Note that direct calculation of the TDE rates from an 
N-body simulation requires appropriate resolution, num¬ 
ber of particles and calibration. We have followed a more 
simplified approach where the N-body simulations only 
provide us with the effective radial density profile at any 
given point in time, from which we then calculate the 
TDE rate semi-analytically, similar to the approach used 
in Mastrobuono-Battisti, Perets & Loeb (2014). This 
avoids the obstacles arising from the direct calculations 
of TDE rates in N-body simulations, at the cost of a less 
accurate result, and the use of the simplifying assump¬ 
tion of an isotropic distribution. 

In the cluster-infall scenario clusters are added to the 
system on equally spaced time intervals. They then in¬ 
spiral and disrupt, and their stars build-up the NSC. To 
a large extent this is captured by the FP models which 
introduce stars in the NSC outskirts; since the infalling 


NSC mass [M©] NSC mass [Mq] 










































Fig. 5. — Evolution of TDE (continuous lines) rates for four different masses of MBHs for models E (left panel) and Sq Sb (right panel) 
for lower rates of SF. The dashed blue lines mark the epoch at which SF stops and the bursts of SF for E and S models respectively. 


X 10 "" 



X lO"'^ 



time [Gyr] 


Fig. 6. — TDEs history for scenarios in which SF that takes place in the inner parsec (0.05-0.1pc), for NSCs that are built up through 
in-situ SF. Left: TDH evaluated for elliptical galaxies (SF occurs only in the first Gyr) for four masses of MBH. Right: TDE histories for 
three different types of galaxies for high SF rates (5 x 10 ~^Mq, similar to Figure[2jl 


clusters disperse in the same distance scales (few pcs), 
they introduce stars in the same manner as assumed in 
these FP models. We therefore compare the result of 
these two scenarios in Figure O As can be seen, both 
scenarios provide very similar results, suggesting that the 
FP calculations can also be used to model the cluster- 
infall scenario. 

4.3. A global model for TDE history 

As described in the previous sections, we have followed 
the evolution of the disruption rates of stars for different 
types of galaxies and evolutionary/SF-history scenarios. 
Combining the rates we obtained for each galaxy type 
with the observational information about the fraction of 
each morphological galaxy type can provide a general 
prediction for the TDH of the universe. This analysis is 
done as follows; for each type of hosted MBH mass we 
obtain the global TDH by integrating the rates we ob¬ 
tained for the different galaxies and MBH mass, weighted 
by their relative fraction 


Fg = / rGi(m_Bg)uGi {rnBH)dMBH (7) 

J Mmin j 

where F Ci is the TDE rate for the i-th galaxy type with 
an MBH of mass msH, and nCiimBn) is its fraction 
of such galaxies in the local universe. The data for the 
fraction of each g alaxy type is taken from observations 
(iCalvi et al.ll201 2n and the M BH mass function is based 
on iKellv and Merlonil (l2012f) . Note that in practice we 
use a sum over the MBH masses for which we have made 
the full calculation, rather than a continuous mass func¬ 
tion. We generally hnd that the TDE rate dependence on 
the MBH mass goes approximately as F^"^^ = 
where 71 = 0.44. We separate the rate dependence re¬ 
sults for each galaxy type, and find that the MBH mass 
goes as Ff^,^ ~ and Ff.^^ ^ for el¬ 

liptical and spiral galaxies respectively. These findings 
of lower TDE rates obtained for higher mass MBHs are 
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x10“® x10'® 




Fig. 7 . — Evolution of TDE rates f or four models with different masses of MBHs, and a SF history which follows the global universal SF 
models IjMadau and Dickinsonlf2Ql‘4l') normalized to the NSCs. (the total number of added stars to each model is equal to the total number 
of stars in models shown in Fig. ^ Left: TDH for NSC evolved from built-up global star formation. Right: TDH for NSC evolved from a 
primordial cusp with global SF. 


X 10“® 



Fig. 8. — Comparison between TDEs history obtained from N- 
body simulations and FP analysis. TDH for an NSC formed 
through clusters infall, compared with the TDE history for the 
FP models with “SF” bursts. The red dashed lines mark the time 
of the infall (N-body) and the blue dashed lines mark the time of 
the SF bursts. 

not surprising and have already been seen and discussed 
in depth in other studies. Previous models had not 
considered continuous SF, however our models for early 
burst of SF (elliptical galaxies) with no later SF are the 
most comparable. Indeed in those cases we find consis¬ 
tent results; ^tde ~ where the range of the 

power law depe ndence found in previous studies is 0.2 < 
71 < O.b le.g. iW ang and Merritt] |2004|: [Merritt||20 13tj : 

IVasiliev and Merrittll201,‘lHStone and Metzgeill2016ll. 

5. SUMMARY 

In our work we provide a first theoretical study of 
the TDE rate history which considers star-formation and 
cluster infall scenarios as models for nuclear stellar clus¬ 
ter formation and evolution. We find that in most scenar¬ 
ios the TDE rates grow with time, due to the growth of 
the NSCs stellar density during their evolution, where the 



Fig. 9. — The integrated cosmological history of TDEs rates per 
galaxy for the built-up and the primordial cusp scenarios for NSC 
formation and evolution. Note that for these final results we in¬ 
tegrated the TDEs rates over MBHs masses between 10® Mq and 
5 X 10’^ Mq. 

TDE rates decrease inversely with the mass of the MBH 
they host. We explored two different origins of the NSC: 
pre-existing NSC cusp that forms very early in the galaxy 
evolution or an NSC that grows more slowly through star 
formation or cluster infalls. NSCs with primordial cusps 
have initially larger stellar populations. The overall TDE 
rates they produce are therefore higher than those calcu¬ 
lated for the cases of NSC which are built-up on longer 
timescales, and they are non-negligible already at early 
times. Consideration of different types of SF histories, 
corresponding to elliptical and spiral galaxies, also show 
a similar trend: the earlier the SF occurs (e.g. in ellipti¬ 
cal where most SF is assumed to occur at the first Gyr), 
the faster the NSC is built, and the faster the TDE rates 
increase. 

We also compared our disruption rates for NSCs with 
SE which occurs in the inner parsec, and found that in 
such SE scenarios a the TDE histories show a different 
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behavior where: during the SF epoch the TDE rates 
rapidly grow, but then decrease and decay to a much 
lower level on a timescale of a few Gyrs. 

Finally, to complement our main results from the 
Fokker-Planck models, we also modeled the TDH using 
data from full N-body simulations of the cluster-infall 
scenario. We found a good match between these rates 
and the FP results corresponding to a source term of 
stars in the NSC outskirts, further confirming the more 
simplified FP analysis. 

Our results are mainly based on simple ID Fokker- 
Planck models. Although these appear to be able to well 
capture the results of N-body simulations (at least for 
the cluster-infall scenario modeled through N-body sim¬ 
ulations), they did not include a detailed galaxy merger 
history, and more complex non-spherical structures of 
NSCs. The role of these aspects and their effects on the 
TDH are not well determined, and should be considered 
in future work 


Future surveys of TDEs (e.g. using the LSST data) 
would be able to provide not only better estimate of the 
TDE rates in the local universe, and their dependence on 
the host galaxy, but also the history of TDEs. As we have 
shown here, the TDH can provide us with information on 
the overall evolution of nuclear stellar clusters and the 
star-formation history and properties in galactic nuclei 
through the combination of the observational data and 
the theoretical studies. 
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FP code used in our simulations. We would also like 
to especially thank the anonymous referee for important 
comments and suggestions that considerably improved 
the manuscript. We acknowledge support from the Tech- 
nion Asher Space Research Institute and the I-CORE 
Program of the Planning and Budgeting Committee and 
The Israel Science Foundation grant 1829/12. 
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